%% bayesian var
function out=bvar_f_w(ps,db)
% modification of bvar_f for matrix operations
% 2020/1/10 
% 


% %% load data
% ps.T=100; % data length
% ps.N=2; % number of variables
% 
% % generate random data for testing
% db.data=randn(ps.N,ps.T); % data: N by T

% ps.prior=2 % minnesota prior
% ps.prior=3 normal-wishart prior
% ps.prior=4 independent wishart prior
%--------------------------------------------------------------------------

%% estimate VAR
% ps.p=2; % number of lags in VAR
% ps.cons=1; % constant in VAR?
% ps.OD=1; % order of MA representation computed
% 
% ps.prior=2; % Minnesota prior
% ps.nsave=10; % draws saved
% ps.nburn=1; % burn in
% ps.ndraws=ps.nsave+ps.nburn;
% ps.nprint=5;
% 
%% make ps.stationary if missing
% disp('ps stationary')
% disp(ps.stationary)

if ~exist('ps.stationary','var')==0
    ps.stationary=0;
end

%% companion form of VAR and OLS estimate 
Y2=db.data';
M=size(Y2,2); % number of variables

% Generate lagged Y matrix. This will be part of the X matrix
Ylag = mlag2(Y2,ps.p); % Y is [T x N]. ylag is [T x (Np)]

% Now define matrix X which has all the R.H.S. variables (constant, lags of
% the dependent variable and exogenous regressors/dummies).

if ps.cons==1
    X1 = [ones(ps.T-ps.p,1) Ylag(ps.p+1:ps.T,:)];
else
    X1 = Ylag(ps.p+1:ps.T,:); 
end
% Get size of final matrix X
[T K] = size(X1);

Y1 = Y2(ps.p+1:ps.T,:); % This is the final Y matrix used for the VAR

Y=Y1;
X=X1;

Z1 = kron(eye(M),X1);
Z=Z1;

%%% OLS estimates 

A_OLS = (X'*X)\(X'*Y); % This is the matrix of regression coefficients
a_OLS = A_OLS(:);         % This is the vector of parameters, i.e. it holds
                          % that a_OLS = vec(A_OLS)

tmp=(Y - X*A_OLS);
                          
SSE = tmp'*tmp;   % Sum of squared errors
SIGMA_OLS = SSE./(T-K+1);


%% initialize

% Initialize Bayesian posterior parameters using OLS values
alpha = a_OLS;     % This is the single draw from the posterior of alpha
ALPHA = A_OLS;     % This is the single draw from the posterior of ALPHA
SSE_Gibbs = SSE;   % This is the SSE based on each draw of ALPHA
SIGMA = SIGMA_OLS; % This is the single draw from the posterior of SIGMA
IXY =  kron(eye(M),(X'*Y));


% Storage space for posterior draws
alpha_draws = zeros(K*M,ps.nsave);   % save draws of alpha
ALPHA_draws = zeros(K,M,ps.nsave);   % save draws of alpha
SIGMA_draws = zeros(M,M,ps.nsave);   % save draws of ALPHA

%% prior hyper parameter 

if ps.prior==2 % Minnesota prior
    A_prior = [zeros(1,M); 1*eye(M); zeros((ps.p-1)*M,M)];  %<---- prior mean of ALPHA (parameter matrix) 
    a_prior = A_prior(:);               %<---- prior mean of alpha (parameter vector)
    
    % Minnesota Variance on VAR regression coefficients
    % First define the hyperparameters 'a_bar_i' 
    % Canova
%     a_bar_1 = 0.2;
%     a_bar_2 = 0.2*0.5;
%     a_bar_3 = 0.2*10^5;
    % Angeletos et al
%     gamma1=0.2;
%     gamma2=0.5;
%     gamma3=2;
%     gamma4=10^5;
    % large BVAR code
%     a_bar_1 = 0.2^2;
%     a_bar_2 = 0.1^2;
%     a_bar_3 = 10^2;
    % original code
    a_bar_1 = 0.2;%0.5; % first lag of its own, smaller = no persistence
    a_bar_2 = 0.1;%0.4; % cross lag of other variables, smaller = no cross variable effect
    a_bar_3 = 0.2*10^4;%10^2; % constant, smaller = tighter

 % ACD 
%     a_bar_1 = 0.1;
%     a_bar_2 = 0.2*0.5;
%     a_bar_3 = 0.2*10^4;
    
    % Now get residual variances of univariate p_MIN-lag autoregressions. Here
    % we just run the AR(p) model on each equation, ignoring the constant
    % and exogenous variables (if they have been specified for the original
    % VAR model)
    p_MIN = 6;
    
% ACD
%      p_MIN = 10;
    
    sigma_sq = zeros(M,1); % vector to store residual variances
    for i = 1:M
        % Create lags of dependent variables   
        Ylag_i = mlag2(Y2(:,i),ps.p);
        Ylag_i = Ylag_i(p_MIN+1:ps.T,:);
        X_i = [ones(ps.T-p_MIN,1) Ylag_i];
        Y_i = Y2(p_MIN+1:ps.T,i);
        % OLS estimates of i-th equation
        alpha_i = (X_i'*X_i)\(X_i'*Y_i);
        sigma_sq(i,1) = (1./(ps.T-p_MIN))*((Y_i - X_i*alpha_i)'*(Y_i - X_i*alpha_i));
    end
    % Now define prior hyperparameters.
    % Create an array of dimensions K x M, which will contain the K diagonal
    % elements of the covariance matrix, in each of the M equations.
    V_i = zeros(K,M);
    
    % index in each equation which are the own lags
    ind = zeros(M,ps.p);
    for i=1:M
        ind(i,:) = ps.cons+i:M:K;
    end
    for i = 1:M  % for each i-th equation
        for j = 1:K   % for each j-th RHS variable
            if ps.cons==1
                if j==1 % if there is constant, use this code
                    V_i(j,i) = a_bar_3*sigma_sq(i,1); % variance on constant
%                     V_i(j,i) = gamma4^2*sigma_sq(i,1); % variance on constant    
                elseif find(j==ind(i,:))>0
                    V_i(j,i) = a_bar_1./(ceil((j-1)/M)^2); % variance on own lags  
%                     V_i(j,i) =(gamma1./(ceil((j-1)/M)^gamma3)); % variance on own lags  
                    % Note: the "ceil((j-1)/M)" command finds the associated lag 
                    % number for each parameter
                else
                    for kj=1:M
                        if find(j==ind(kj,:))>0
                            ll = kj;                   
                        end
                    end                 % variance on other lags  
                    V_i(j,i) = (a_bar_2*sigma_sq(i,1))./((ceil((j-1)/M)^2)*sigma_sq(ll,1)); 
%                     V_i(j,i) = (gamma1*gamma2*sigma_sq(i,1)./((ceil((j-1)/M)^gamma3)*sigma_sq(ll,1)));           
                end
            else   % if no constant is defined, then use this code
                if find(j==ind(i,:))>0
                    V_i(j,i) = a_bar_1./(ceil(j/M)^2); % variance on own lags
%                     V_i(j,i) =(gamma1./(ceil((j-1)/M)^gamma3)); % variance on own lags  

                else
                    for kj=1:M
                        if find(j==ind(kj,:))>0
                            ll = kj;
                        end                        
                    end                 % variance on other lags  
                    V_i(j,i) = (a_bar_2*sigma_sq(i,1))./((ceil(j/M)^2)*sigma_sq(ll,1));            
%                     V_i(j,i) = (gamma1*gamma2*sigma_sq(i,1)./((ceil((j-1)/M)^gamma3)*sigma_sq(ll,1)));           
                end
            end
        end
    end
    
    % Now V is a diagonal matrix with diagonal elements the V_i
    V_prior = diag(V_i(:));  % this is the prior variance of the vector alpha
    
    %NOTE: No prior for SIGMA. SIGMA is simply a diagonal matrix with each
    %diagonal element equal to sigma_sq(i). See Kadiyala and Karlsson (1997)
    SIGMA = diag(sigma_sq);
    
elseif ps.prior==2.5 % ACD Minnesota
    A_prior = [zeros(1,M); 1*eye(M); zeros((ps.p-1)*M,M)];  %<---- prior mean of ALPHA (parameter matrix) 
    a_prior = A_prior(:);               %<---- prior mean of alpha (parameter vector)
   
    lb1=0.2;
    lb2=0.5;
%     lb3=2;
%     lb4=10^5;
%     rho=ones(M);
%     se=zeros(M);
        sigma_sq = zeros(M,1); % vector to store residual variances
        sigma_sq1 = zeros(M,1); % vector to store residual variances
    for i = 1:M
        % Create lags of dependent variables   
        Ylag_i = Y2(1:ps.T-1,i); %mlag2(Y2(:,i),ps.p);
        %Ylag_i = Ylag_i(p_MIN+1:ps.T,:);
        X_i = [ones(ps.T-1,1) Ylag_i];
        Y_i = Y2(2:ps.T,i);
        % OLS estimates of i-th equation
        alpha_i = (X_i'*X_i)\(X_i'*Y_i);
        sigma_sq(i,1) = (1./(size(Y_i,1)-2))*((Y_i - X_i*alpha_i)'*(Y_i - X_i*alpha_i));
        sigma_sq1(i,1) = sqrt(sigma_sq(i,1));
    end
        V_i = zeros(K,M);
    
%      H=zeros(K,M);
%      TH=zeros(K,M);
%      for i=1:M
%          k=1;
%          for l=1:ps.p
%              for j=1:M
%                  if j==i
%                      tmp=lb1/l^lb3;
%                      H(k,i)=tmp*tmp;
%                      if l==1
%                          TH(k,i)=rho(i,1);
%                      end
%                  else
%                      tmp=sigma_sq1(i,1)*lb1*lb2/(sigma_sq1(j,1)*l^lb3);
%                      H(k,i)=tmp*tmp;
%                  end
%                  k=k+1;
%              end
%          end
%          tmp=sigma_sq1(i,1)*lb4;
%          H(k,i)=tmp*tmp;
%          TH(k,i)=0;
%      end
    a_bar_1 = 0.7;%0.5; % first lag of its own, smaller = no persistence
    a_bar_2 = 0.4; %lb1*lb2;%0.4; % cross lag of other variables, smaller = no cross variable effect
    a_bar_3 = 10^3;%10^5;%10^2; % constant, smaller = tighter
    
    % index in each equation which are the own lags
    ind = zeros(M,ps.p);
    for i=1:M
        ind(i,:) = ps.cons+i:M:K;
    end
    for i = 1:M  % for each i-th equation
        for j = 1:K   % for each j-th RHS variable
            if ps.cons==1
                if j==1 % if there is constant, use this code
                    V_i(j,i) = (a_bar_3*sigma_sq1(i,1))^2; % variance on constant
%                     V_i(j,i) = gamma4^2*sigma_sq(i,1); % variance on constant    
                elseif find(j==ind(i,:))>0
                    V_i(j,i) = (a_bar_1./(ceil((j-1)/M)^2))^2; % variance on own lags  
%                     V_i(j,i) =(gamma1./(ceil((j-1)/M)^gamma3)); % variance on own lags  
                    % Note: the "ceil((j-1)/M)" command finds the associated lag 
                    % number for each parameter
                else
                    for kj=1:M
                        if find(j==ind(kj,:))>0
                            ll = kj;                   
                        end
                    end                 % variance on other lags  
                    V_i(j,i) = ((a_bar_2*sigma_sq1(i,1))./((ceil((j-1)/M)^2)*sigma_sq1(ll,1)))^2; 
%                     V_i(j,i) = (gamma1*gamma2*sigma_sq(i,1)./((ceil((j-1)/M)^gamma3)*sigma_sq(ll,1)));           
                end
            else   % if no constant is defined, then use this code
                if find(j==ind(i,:))>0
                    V_i(j,i) = a_bar_1./(ceil(j/M)^2); % variance on own lags
%                     V_i(j,i) =(gamma1./(ceil((j-1)/M)^gamma3)); % variance on own lags  

                else
                    for kj=1:M
                        if find(j==ind(kj,:))>0
                            ll = kj;
                        end                        
                    end                 % variance on other lags  
                    V_i(j,i) = (a_bar_2*sigma_sq(i,1))./((ceil(j/M)^2)*sigma_sq(ll,1));            
%                     V_i(j,i) = (gamma1*gamma2*sigma_sq(i,1)./((ceil((j-1)/M)^gamma3)*sigma_sq(ll,1)));           
                end
            end
        end
    end
            
        
    
    
    % Now V is a diagonal matrix with diagonal elements the V_i
    V_prior = diag(V_i(:));  % this is the prior variance of the vector alpha
    %NOTE: No prior for SIGMA. SIGMA is simply a diagonal matrix with each
    %diagonal element equal to sigma_sq(i). See Kadiyala and Karlsson (1997)
    % SIGMA = diag(sigma_sq);
    SIGMA= SIGMA_OLS; %SSE./(T); %SIGMA_OLS;   
    
elseif ps.prior==3 % Normal-Whishart
    % Hyperparameters on a ~ N(a_prior, SIGMA x V_prior) 
    A_prior = 0*ones(K,M);   %<---- prior mean of ALPHA (parameter matrix)
    a_prior = A_prior(:);    %<---- prior mean of alpha (parameter vector)
    V_prior = 1*eye(K);%10*eye(K);     %<---- prior variance of alpha
    
    % Hyperparameters on inv(SIGMA) ~ W(v_prior,inv(S_prior))
    v_prior = M+1;                 %<---- prior Degrees of Freedom (DoF) of SIGMA
    S_prior = 1*eye(M);%1*eye(M);      %<---- prior scale of SIGMA
    inv_S_prior=eye(length(S_prior))/S_prior;
    
 elseif ps.prior == 4  % Independent Normal-Wishart
    n = K*M; % Total number of parameters (size of vector alpha)
  %  A_prior = [zeros(1,M); 1*eye(M); zeros((ps.p-1)*M,M)];  %<---- prior mean of ALPHA (parameter matrix)
    a_prior = 0*ones(n,1);   %<---- prior mean of alpha (parameter vector)
   % a_prior = A_prior(:);    %<---- prior mean of alpha (parameter vector)
    V_prior = 10*eye(n);     %<---- prior variance of alpha
    
    % Hyperparameters on inv(SIGMA) ~ W(v_prior,inv(S_prior))
    v_prior = M+1;             %<---- prior Degrees of Freedom (DoF) of SIGMA
    S_prior = eye(M);            %<---- prior scale of SIGMA
    inv_S_prior = inv(S_prior);

end

%% start sampling
%tic
%disp('Number of iterations');
for irep = 1:ps.ndraws  %Start the Gibbs "loop"
   
    %for irep = 1:ps.ndraws  %Start the Gibbs "loop"
 %   if mod(irep,ps.nprint) == 0 % print iterations
%         disp(irep);
 %       toc;
 %   end

    if ps.prior==2 
    %--------- Draw ALPHA and SIGMA with Minnesota Prior
        for i = 1:M
 %           V_post = inv( inv(V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K)) + inv(SIGMA(i,i))*(X'*X) );
            V_post_tmp=inv(V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K)) + (SIGMA(i,i))\(X'*X);
            V_post = eye(length(V_post_tmp))/V_post_tmp;
            
%            a_post = V_post*(inv(V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K))*a_prior((i-1)*K+1:i*K,1) + inv(SIGMA(i,i))*X'*Y(:,i));
            a_post = V_post*((V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K))\a_prior((i-1)*K+1:i*K,1) + (SIGMA(i,i))\(X'*Y(:,i)));

            alpha((i-1)*K+1:i*K,1) = a_post + chol(V_post)'*randn(K,1); % Draw alpha
        end
        ALPHA = reshape(alpha,K,M); % Create draw in terms of ALPHA
        
        % SIGMA in this case is a known matrix, whose form is decided in
        % the prior (see prior_hyper.m)
    %--------- Draw ALPHA and SIGMA with Normal-Wishart Prior
    elseif ps.prior==2.5 
    %--------- Draw ALPHA and SIGMA with Minnesota Prior    
    
    %struct obj_minnesota_priors
%  	th0 coefficient
% 	sigma0 V_prior
% 	isigma0 V_prior\eye
% 	gamma0 eye(M)
% 	t0 M+1
            %DRAW COEFFICIENTS
            ISIG0=V_prior\eye(length(V_prior));
%            ISIG0=(V_prior+1e-7*eye(length(V_prior)))\eye(length(V_prior));
            
            isig = inv(chol(SIGMA'));
            ISIG = isig*isig';
            tmp = (ISIG0 + kron(ISIG,X'*X));
            tmp = (tmp+tmp')/2;
            itmp = inv(chol(tmp));
            Vstar = itmp*itmp';
            Vstar = (Vstar+Vstar')/2;
            Mstar   = Vstar*(ISIG0*a_prior+kron(ISIG,(X'*X))*A_OLS(:));
            chk=0;
            while chk==0
                alpha=Mstar+chol(Vstar)'*randn(M*(ps.p*M+1),1);
                ALPHA=reshape(alpha,K,M);
            %companion form coefficient hx : X(t)=hx*X(t-1)+netashock*e(t)
            hx=zeros(ps.N*ps.p,ps.N*ps.p); % preallocation
            hx(ps.N+1:end,1:ps.N*(ps.p-1))=eye((ps.p-1)*ps.N);

            tmp=ALPHA(2:end,:); % take out constants
            hx(1:ps.N,1:ps.N*ps.p)=transpose(tmp); % transpose 
            m=max(abs(eig(hx)));  
%             disp(m)
             if m<0.999999
                chk=1;
             end
            end
            
            % draw Volatility conditional on new coefficients
%             resid=(Y-X*ALPHA)-mean(Y-X*ALPHA,1).*ones(T,1);
%             T1=M+1+T;
%             GAMMA1=eye(M)+resid'*resid;
%             IG1=GAMMA1\eye(M);
%             IG1=(IG1+IG1')/2;
%             z=randn(T1,M)*chol(IG1);
%             SIGMA=(z'*z)\eye(M);
            
         %    disp(irep)
            
%        if ps.stationary==1
%         % Posterior of alpha|SIGMA,Data ~ Normal
%          chk=0;
%         while chk==0
%             for i = 1:M
%             V_post_tmp=inv(V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K)) + (SIGMA(i,i))\(X'*X);
%             V_post = eye(length(V_post_tmp))/V_post_tmp;
%             
% %            a_post = V_post*(inv(V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K))*a_prior((i-1)*K+1:i*K,1) + inv(SIGMA(i,i))*X'*Y(:,i));
%             a_post = V_post*((V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K))\a_prior((i-1)*K+1:i*K,1) + (SIGMA(i,i))\(X'*Y(:,i)));
% 
%             alpha((i-1)*K+1:i*K,1) = a_post + chol(V_post)'*randn(K,1); % Draw alpha
%             end
%             ALPHA = reshape(alpha,K,M); % Create draw in terms of ALPHA
%           % companion form coefficient hx : X(t)=hx*X(t-1)+netashock*e(t)
%             hx=zeros(ps.N*ps.p,ps.N*ps.p); % preallocation
%             hx(ps.N+1:end,1:ps.N*(ps.p-1))=eye((ps.p-1)*ps.N);
% 
%             tmp=ALPHA(2:end,:); % take out constants
%             hx(1:ps.N,1:ps.N*ps.p)=transpose(tmp); % transpose 
%             m=max(abs(eig(hx)));
%            % disp(m)
%             if m<0.999999
%                 chk=1;
%             end
%         end
%      elseif ps.stationary==0
%         for i = 1:M
%  %           V_post = inv( inv(V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K)) + inv(SIGMA(i,i))*(X'*X) );
%             V_post_tmp=inv(V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K)) + (SIGMA(i,i))\(X'*X);
%             V_post = eye(length(V_post_tmp))/V_post_tmp;
%             
% %            a_post = V_post*(inv(V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K))*a_prior((i-1)*K+1:i*K,1) + inv(SIGMA(i,i))*X'*Y(:,i));
%             a_post = V_post*((V_prior((i-1)*K+1:i*K,(i-1)*K+1:i*K))\a_prior((i-1)*K+1:i*K,1) + (SIGMA(i,i))\(X'*Y(:,i)));
% 
%             alpha((i-1)*K+1:i*K,1) = a_post + chol(V_post)'*randn(K,1); % Draw alpha
%         end
%         ALPHA = reshape(alpha,K,M); % Create draw in terms of ALPHA
% 
%       end

    elseif ps.prior == 3
        % ******Get all the required quantities for the posteriors   
        inv_V_prior=eye(length(V_prior))/V_prior;
        V_post_tmp=inv_V_prior + (X'*X);
        V_post = eye(length(V_post_tmp))/V_post_tmp;
        %A_post = V_post*(inv(V_prior)*A_prior + (X'*X)*A_OLS);
        A_post = V_post*((V_prior)\A_prior + (X'*X)*A_OLS);
        
        a_post = A_post(:);
   
        %S_post = SSE + S_prior + A_OLS'*(X'*X)*A_OLS + A_prior'*inv(V_prior)*A_prior - A_post'*(inv(V_prior) + (X'*X))*A_post;
     S_post = SSE + S_prior + A_OLS'*(X'*X)*A_OLS + A_prior'*inv_V_prior*A_prior - A_post'*(inv_V_prior + (X'*X))*A_post;

        v_post = T + v_prior;
    
        % This is the covariance for the posterior density of alpha
        COV = kron(SIGMA,V_post);
        
        if ps.stationary==1
        % Posterior of alpha|SIGMA,Data ~ Normal
         chk=0;
        while chk==0
         alpha = a_post + chol(COV)'*randn(K*M,1);  % Draw alpha
         ALPHA = reshape(alpha,K,M); % Draw of ALPHA
          % companion form coefficient hx : X(t)=hx*X(t-1)+netashock*e(t)
            hx=zeros(ps.N*ps.p,ps.N*ps.p); % preallocation
            hx(ps.N+1:end,1:ps.N*(ps.p-1))=eye((ps.p-1)*ps.N);

            tmp=ALPHA(2:end,:); % take out constants
            hx(1:ps.N,1:ps.N*ps.p)=transpose(tmp); % transpose 
            m=max(abs(eig(hx)));
           % disp(m)
            if m<0.999999
                chk=1;
            end
        end
        elseif ps.stationary==0
          alpha = a_post + chol(COV)'*randn(K*M,1);  % Draw alpha
          ALPHA = reshape(alpha,K,M); % Draw of ALPHA

        end
        
        % Posterior of SIGMA|ALPHA,Data ~ iW(inv(S_post),v_post)
        inv_S_post=eye(length(S_post))/S_post;
        wish_tmp=wish(inv_S_post,v_post);
        SIGMA = eye(length(wish_tmp))/wish_tmp;% Draw SIGMA
        
     elseif ps.prior == 4
       %  disp(irep)
          chk=0;
        while chk==0
        VARIANCE = kron(inv(SIGMA),eye(T));
        V_post_tmp = V_prior + Z'*VARIANCE*Z;
        V_post = eye(length(V_post_tmp))/V_post_tmp;
        %V_post = inv(V_prior + Z'*VARIANCE*Z);
        a_post = V_post*(V_prior*a_prior + Z'*VARIANCE*Y(:));
        alpha = a_post + chol(V_post)'*randn(n,1); % Draw of alpha
        
        ALPHA = reshape(alpha,K,M); % Draw of ALPHA
        % make sure it's stationary
            hx=zeros(ps.N*ps.p,ps.N*ps.p); % preallocation
            hx(ps.N+1:end,1:ps.N*(ps.p-1))=eye((ps.p-1)*ps.N);

            tmp=ALPHA(2:end,:); % take out constants
            hx(1:ps.N,1:ps.N*ps.p)=transpose(tmp); % transpose 
            m=max(abs(eig(hx)));
           % disp(m)
            if m<0.999999
                chk=1;
            end
        end
        % Posterior of SIGMA|ALPHA,Data ~ iW(inv(S_post),v_post)
        v_post = T + v_prior;
        S_post = S_prior + (Y - X*ALPHA)'*(Y - X*ALPHA);
%        SIGMA = inv(wish(inv(S_post),v_post));% Draw SIGMA
        
        S_post_tmp = eye(length(S_post))/S_post;
        wish_tmp = wish(S_post_tmp,v_post);
        SIGMA = eye(length(wish_tmp))/wish_tmp;% Draw SIGMA
  
    end
        if irep > ps.nburn
        
        alpha_draws(:,irep-ps.nburn) = alpha;
        ALPHA_draws(:,:,irep-ps.nburn) = ALPHA;
        SIGMA_draws(:,:,irep-ps.nburn) = SIGMA;

        
        end


    
end

%% output


out.ALPHA_draws=ALPHA_draws;
out.SIGMA_draws=SIGMA_draws;

out.ALPHA_OLS=A_OLS;
out.SIGMA_OLS=SIGMA_OLS;

